rm(list = ls())
library(magrittr)
# 生成模拟数据
set.seed(1024)
y1 <- rnorm(70, mean = 0.5, sd = 0.3)
y2 <- rnorm(30, mean = 2, sd = 0.6)
y <- c(y1,y2)

# 估计
# 初始化
pi <- 0.5
theta <- c(u1 = y1[1], sgm1 = sd(y), u2 = y1[2], sgm2 = sd(y))
gamma <- pi * dnorm(y, mean = theta[3], sd = theta[4])/
  ((1-pi) * dnorm(y, mean = theta[1], sd = theta[2]) + pi * dnorm(y, mean = theta[3], sd = theta[4]))

# 似然函数
lhfun <- function(theta, ydata, gamma,pi){
  likelihood <- (1-gamma)*log(dnorm(ydata, mean = theta[1],sd = theta[2]))+
    gamma*log(dnorm(ydata, mean = theta[3], sd = theta[4])) +
    (1-gamma) * log(1-pi) + gamma*log(pi)
  return(-sum(likelihood))
}

# 迭代结果存在estpara中
estpara <- NULL
eps <- 0.5
while (any(eps > 0.0001)) {
  # E步
  if (!is.null(estpara)) gamma <- pi * dnorm(y, mean = estpara[nrow(estpara),3], sd = estpara[nrow(estpara),4])/
    ((1-pi) * dnorm(y, mean = estpara[nrow(estpara),1], sd = estpara[nrow(estpara),2]) + 
  pi * dnorm(y, mean = estpara[nrow(estpara),3], sd = estpara[nrow(estpara),4]))
  pi <- mean(gamma)
  # M 步
  medpar <- optimx::optimx(par = theta, fn = lhfun, ydata = y, gamma = gamma,pi=pi, method = 'BFGS') %>% 
    .[1:4] %>% as.data.frame()
  estpara <- rbind(estpara,medpar)
  # 最近两次迭代的差距
  if (nrow(estpara) >= 2) eps <- abs(estpara[nrow(estpara),] - estpara[nrow(estpara)-1,])
}
